home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Language/OS - Multiplatform Resource Library
/
LANGUAGE OS.iso
/
cpp_libs
/
newmat03.lha
/
newmat03
/
newmat1.cxx
< prev
next >
Wrap
C/C++ Source or Header
|
1993-08-08
|
8KB
|
283 lines
//$$ newmat1.cxx Utilities and matrix type functions
// Copyright (C) 1991: R B Davies and DSIR
#define WANT_STREAM // to include stream package
#include "include.hxx"
#include "newmat.hxx"
/**************************** error handlers *******************************/
void MatrixError(char* message)
{
cerr << "\nError detected by matrix package: " << message << "\n";
exit(0);
}
void MatrixErrorNoSpace(void* v)
// give error message if v is null
{
if (v) return;
cerr << "\nError detected by matrix package: no space\n";
exit(0);
}
/************************* ExeCounter functions *****************************/
// the next statment may need to be deleted form some compilers
int ExeCounter::nreports = 0;
ExeCounter::ExeCounter(int xl, int xf) : line(xl), fileid(xf), nexe(0) {}
ExeCounter::~ExeCounter()
{
nreports++;
cout << nreports << " " << fileid << " " << line << " " << nexe << "\n";
}
/************************* MatrixType functions *****************************/
BOOL MatrixType::operator==(const MatrixType& t) const
{ return (type == t.type); }
BOOL MatrixType::operator!=(const MatrixType& t) const
{ return (type != t.type); }
MatrixType MatrixType::operator+(const MatrixType& mt) const
{
Type type1=mt.type;
if (type==UnSp || type1==UnSp) return UnSp;
if (type==type1) return type; // won't work for orthogonal
if (type==ColV || type1==ColV) return ColV;
if (type==RowV || type1==RowV) return RowV;
if (type==EqEl && (type1==Sym || type1==Diag)) return Sym;
if (type1==EqEl && (type==Sym || type==Diag)) return Sym;
if (type==EqEl || type1==EqEl) return Rect;
if (type==Diag) return type1; // won't work for assym
if (type1==Diag) return type;
return Rect;
}
MatrixType MatrixType::operator-(const MatrixType& mt) const
{
return *this+mt; // won't work for orthogonal or pos def
}
MatrixType MatrixType::operator*(const MatrixType& mt) const
{
Type type1=mt.type;
if (type==UnSp || type1==UnSp) return UnSp;
if (type1==ColV) { if (type==RowV) return Diag; else return ColV; }
if (type==RowV) return RowV;
if (type==Sym || type1==Sym) return Rect;
if (type==type1) return type; // won't work for pos def or assym
if (type==EqEl || type1==EqEl) return Rect;
if (type==Diag) return type1;
if (type1==Diag) return type;
return Rect;
}
BOOL MatrixType::operator>=(const MatrixType& mt) const
{
Type type1=mt.type;
if (type==type1) return TRUE; // by definition
if (type==Rect || type==RowV || type==ColV) return TRUE;
if (type==EqEl) return FALSE;
if (type1==Diag) return TRUE; // won't work for pos def or assym
return FALSE;
}
MatrixType MatrixType::i() const
{
switch (type)
{
case UnSp: return UnSp;
case UT: return UT;
case LT: return LT;
case Rect: return Rect;
case Sym: return Sym;
case Diag: return Diag;
case RowV: return Diag;
case ColV: return Diag;
case EqEl: return Diag;
case Crout: return UnSp;
}
}
MatrixType MatrixType::t() const
{
switch (type)
{
case UnSp: return UnSp;
case UT: return LT;
case LT: return UT;
case Rect: return Rect;
case Sym: return Sym;
case Diag: return Diag;
case RowV: return ColV;
case ColV: return RowV;
case EqEl: return EqEl;
case Crout: return UnSp;
}
}
MatrixType MatrixType::sub() const
{
switch (type)
{
case UnSp: return UnSp;
case UT: return Rect;
case LT: return Rect;
case Rect: return Rect;
case Sym: return Rect;
case Diag: return Rect;
case RowV: return RowV;
case ColV: return ColV;
case EqEl: return EqEl;
case Crout: return UnSp;
}
}
MatrixType MatrixType::ssub() const
{
return *this; // won't work for selection matrix
}
MatrixType::operator char*() const
{
switch (type)
{
case UnSp: return "UnSp ";
case UT: return "UT ";
case LT: return "LT ";
case Rect: return "Rect ";
case Sym: return "Sym ";
case Diag: return "Diag ";
case RowV: return "RowV ";
case ColV: return "ColV ";
case EqEl: return "EqEl ";
case Crout: return "Crout";
}
}
GeneralMatrix* MatrixType::New() const
{
GeneralMatrix* gm;
switch (type)
{
case UnSp: MatrixError("Can't build matrix type UnSp");
case UT: gm = new UpperTriangularMatrix; break;
case LT: gm = new LowerTriangularMatrix; break;
case Rect: gm = new Matrix; break;
case Sym: gm = new SymmetricMatrix; break;
case Diag: gm = new DiagonalMatrix; break;
case RowV: gm = new RowVector; break;
case ColV: gm = new ColumnVector; break;
case EqEl: MatrixError("Can't build matrix type EqEl");
case Crout: MatrixError("Can't build matrix type Crout");
}
MatrixErrorNoSpace(gm);
gm->Protect(); return gm;
}
GeneralMatrix* MatrixType::New(int nr, int nc) const
{
GeneralMatrix* gm;
switch (type)
{
case UnSp: MatrixError("Can't build matrix type UnSp");
case UT: gm = new UpperTriangularMatrix(nr); break;
case LT: gm = new LowerTriangularMatrix(nr); break;
case Rect: gm = new Matrix(nr, nc); break;
case Sym: gm = new SymmetricMatrix(nr); break;
case Diag: gm = new DiagonalMatrix(nr); break;
case RowV: if (nr!=1) MatrixError("Illegal Row Vector");
gm = new RowVector(nc); break;
case ColV: if (nc!=1) MatrixError("Illegal Column Vector");
gm = new ColumnVector(nr); break;
case EqEl: MatrixError("Can't build matrix type EqEl");
case Crout: MatrixError("Can't build matrix type Crout");
}
MatrixErrorNoSpace(gm);
gm->Protect(); return gm;
}
void TestTypeAdd()
{
MatrixType list[] = { MatrixType::UT,
MatrixType::LT,
MatrixType::Rect,
MatrixType::Sym,
MatrixType::Diag,
MatrixType::RowV,
MatrixType::ColV,
MatrixType::EqEl };
cout << "+ ";
for (int i=0; i<MatrixType::nTypes(); i++) cout << (char*)list[i] << " ";
cout << "\n";
for (i=0; i<MatrixType::nTypes(); i++)
{
cout << (char*)(list[i]) << " ";
for (int j=0; j<MatrixType::nTypes(); j++)
cout << (char*)(list[j]+list[i]) << " ";
cout << "\n";
}
cout << "\n";
}
void TestTypeMult()
{
MatrixType list[] = { MatrixType::UT,
MatrixType::LT,
MatrixType::Rect,
MatrixType::Sym,
MatrixType::Diag,
MatrixType::RowV,
MatrixType::ColV,
MatrixType::EqEl };
cout << "* ";
for (int i=0; i<MatrixType::nTypes(); i++)
cout << (char*)list[i] << " ";
cout << "\n";
for (i=0; i<MatrixType::nTypes(); i++)
{
cout << (char*)list[i] << " ";
for (int j=0; j<MatrixType::nTypes(); j++)
cout << (char*)(list[j]*list[i]) << " ";
cout << "\n";
}
cout << "\n";
}
void TestTypeOrder()
{
MatrixType list[] = { MatrixType::UT,
MatrixType::LT,
MatrixType::Rect,
MatrixType::Sym,
MatrixType::Diag,
MatrixType::RowV,
MatrixType::ColV,
MatrixType::EqEl };
cout << ">= ";
for (int i = 0; i<MatrixType::nTypes(); i++)
cout << (char*)list[i] << " ";
cout << "\n";
for (i=0; i<MatrixType::nTypes(); i++)
{
cout << (char*)list[i] << " ";
for (int j=0; j<MatrixType::nTypes(); j++)
cout << ((list[j]>=list[i]) ? "Yes " : "No ");
cout << "\n";
}
cout << "\n";
}